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We propose a unified physical framework for transport in variably saturated porous media. This 
approach allows fluid flow and solute migration to be treated as ensemble averages of fluid and 
solute particles, respectively. We consider the cases of homogeneous and heterogeneous porous 
materials. Within a fractal mobile-immobile (MIM) continuous time random walk framework, the 
heterogeneity will be characterized by algebraically decaying particle retention-times. We derive 
the corresponding (nonlinear) continuum limit partial differential equations and we compare their 
solutions to Monte Carlo simulation results. The proposed methodology is fairly general and can be 
used to track fluid and solutes particles trajectories, for a variety of initial and boundary conditions. 



I. INTRODUCTION 

Accurately predicting the spreading of a chemical 
species through unsaturated porous materials (i.e., ma- 
terials with locally-varying fluid content) is key to mas- 
tering such technological challenges as polluted sites re- 
mediation, environmental protection, and waste manage- 
ment P, H, H, 0, [E ■ The evolution of the solute concen- 
tration profile in the traversed medium cannot be a priori 
decoupled from that of the fluid flow, which is ultimately 
responsible for the advection and dispersion mechanisms 
of the solutes. Moreover, the effects of the fluid distri- 
bution are often combined with those of spatial hetero- 
geneities. Such heterogeneities can be present both at the 
pore-scale (microscopic) and at the Darcy-scale (macro- 
scopic). As a consequence, experimental results reported 
in literature often show that solutes concentration dis- 
plays non-Fickian features, such as breakthrough curves 
with long tails and non-Gaussian spatial profiles @, 0, § • 

In this respect, there exists an increasing need for re- 
liable numerical techniques to tackle flow and transport 
problems. In this work, we address the issue of determin- 
ing the fluid content and the solute concentration pro- 
files within unsaturated homogeneous as well as hetero- 
geneous media by resorting to a random walk approach. 
Random walks are extensively used to describe solutes 
transport in saturated media (j| [l^ , although their ap- 
plication to unsaturated flows appears to be somehow 
neglected While complementing each other, the 

random walk and the continuum-limit approaches dis- 
play specific advantages and disadvantages. Random 
walks, for instance, do not introduce the spurious nu- 
merical dispersion typical of Eulerian (continuum) nu- 
merical schemes 0, [HI [HI • As such, random walks are 
particularly well suited to deal with unsaturated mate- 
rials, where sharp contrasts between stagnant and fluid- 
saturated regions, or at macroscopic heterogeneity inter- 
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faces, may give rise to steep propagating fronts. Eulerian 
(continuum) numerical schemes, on the other hand, are 
generally faster than the corresponding Monte Carlo sim- 
ulations. 

We begin our analysis by illustrating flow and trans- 
port in homogeneous media, and detail how the Richards 
equation for the fluid flow and the Advection-Dispersion 
Equation (ADE) for the solutes can be recast in a Fokkcr- 
Planck Equation (FPE) form. The FPE governs the 
probability density (pdf) of finding a walker (a fluid or 
solute parcel, respectively) at a given position at a given 
time. The central idea is that these walkers perform 
stochastic trajectories in the traversed medium. Taking 
the ensemble average of the fluid and solute parcels tra- 
jectories yields the desired macroscopic quantities, i.e., 
the fluid content and solutes concentration profiles. As 
the fluid movement and the solute transport are (nonlin- 
ear ly) coupled via the macroscopic governing equations, 
the underlying stochastic trajectories also display a non- 
linear coupling. On the other hand, the homogeneity hy- 
pothesis ensures that the trajectories carry no memory 
of the past, so that the particles dynamics is Markovian 

Then, we focus our attention on unsaturated hetero- 
geneous materials, where non-Fickian behaviors are en- 
hanced by the interplay between nonlincarities in the 
flow patterns and complex spatial structures: the rel- 
ative strength of these processes determines the precise 
details of the solutes distribution. We model the effects of 
complex nonhomogeneous spatial structures by introduc- 
ing the possibility of trapping events between successive 
particles displacements, as customary within a continu- 
ous time random walk (CTRW) approach [H, [H|]. The 
resulting broad distribution of waiting times at each vis- 
ited site characterizes the broad velocity spectrum that 
is often observed in heterogeneous media. 

This paper is organized as follows: in Sec. [TTJ, we re- 
vise the physical equations that govern the coupled flow- 
transport problem for nonstationary variable-saturation 
conditions. Then, in Sec. IIIII we present a general non- 
linear random walk approach to the simulation of fluid 
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FIG. 1: Saturation level spatial profiles s(x,t) at fixed time 
t = 1, for varying \. Numerical integration of Eq. (j4]) is 
displayed as solid lines, Monte Carlo simulation as symbols: 
X = 0.1 (dots), x = 0.5 (circles), x = 1 (squares), and x = 2 
(crosses). The other coefficients are: Ko = 0.05, vo = 0.3, and 
v = 0.4. 



FIG. 2: Saturation level spatial profiles s(x,t) at fixed time 
t — 1, for varying ia Numerical integration of Eq. (j4| is 
displayed as solid lines, Monte Carlo simulation as symbols: 
v — (dots), v — 0.5 (circles), v — 1 (squares), and i/ = 2 
(crosses). The other coefficients are: ko = 0.05, vq = 0.3, and 
X = 0.5. 



and contaminant particles in locally homogeneous me- 
dia. These simulation schemes are compared with nu- 
merical solutions of the governing equations in Sec. ITVl 
for a variety of initial and boundary conditions. The 
case of discontinuous physical properties of the traversed 
media is addressed, as well. In Sec. |V] we extend our 
results to heterogeneous porous media, and derive the 
corresponding macroscopic governing equations. Finally, 
conclusions are drawn in Sec. IVII 



II. GOVERNING EQUATIONS IN 
UNSATURATED HOMOGENEOUS MEDIA 

In the following, we shall briefly recall the physical 
equations that govern non-stationary flow-transport pro- 
cesses in unsaturated homogeneous porous media. Con- 
sider a vertical column of length I and radius r <C £, 
so that the flow-transport process can be considered one- 
dimensional along the longitudinal direction. Let the col- 
umn be filled with a homogeneous porous material, and 
suppose that the medium has an initial variable satura- 
tion. The fluid flow dynamics within such region can be 
described in terms of the (dimensionless) volumetric fluid 
content < s(x,t) < 1 [!, [Hj], whose evolution is ruled 
by the continuity equation 

d t s(x,t) = -d x j s (x,t), (1) 

provided that the porosity is constant, i.e., the soil 
skeleton is rigid [Ta]. When s(x,t) = 1 everywhere, 
the medium is fully saturated in fluid. The so-called 
Buckingham-Darcy flux (or generalized Darcy's law) 



j s (x,t) [L/T] is provided by the constitutive equation 

j s (x,t) = K(s)[l-d x h(s)}, (2) 

where the quantity K(s) [L/T] is the saturation- 
dependent hydraulic conductivity, and h(s) [L] is the 
saturation-dependent capillary pressure (3l.ll2j. Then, by 
introducing the capillary diffusivity k(s) = K(s)dh/ds 
[L 2 /T], we can combine Eqs. ([J) and ((2|) so to obtain 

d t s(x, t) = d x K(s)d x s(x, t) - d x K(s). (3) 

In this form, Eq. Q has been first derived by 
Richards [l6|. The Richards equation has been exten- 
sively adopted in describing the dynamics of fluid flows 
during wetting processes in soils (see the discussion in [l5j 
and references therein). Without loss of generality, we 
can finally put this equation in conservative form by con- 
veniently defining a 'velocity' v(s) = K(s)/s [L/T]: 

d t s(x,t) = -d x [v(s) - K(s)d x ]s(x,t). (4) 

The term v(s) plays the role of an effective velocity for 
the fluid particles and represents the gravitational contri- 
bution to the flow dynamics. Remark that equation 
is nonlinear, in that the diffusion n(s) and advection 
v(s) coefficients depend in general on s. In some spe- 
cial cases, it is nonetheless possible to obtain analytical 
solutions, b y re sorting to the scaled (Boltzmann) variable 
e^xt- 1 / 2 [Hp. 

Flow dynamics must be supplemented by the initial 
and boundary conditions. To set the ideas, as a repre- 
sentative example we may impose s(0, t) — 1, i.e., we 
keep the inlet on the column at a constant full satu- 
ration. This condition may be physically achieved by 
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FIG. 3: Breakthrough curves j s (£,t) as a function of time 
t for discontinuous k coefficient. The interface is located at 
Xd — £/2. Numerical integration is displayed as solid lines, 
Monte Carlo simulation as symbols: kq — 0.4 in the left layer 
and Ko = 0.05 in the right layer (crosses); ko = 0.05 in the 
right layer and ko — 0.4 in the left layer (circles). The other 
coefficients are: v — 0.2, vq = 0.3, and \ — 0.1. 

putting the porous column in contact with a fluid reser- 
voir (infiltration process). Along the column, we initially 
assign a given saturation distribution, for instance a con- 
stant profile s(x,0) = so, for x > 0. Finally, at the 
outlet of the column we prescribe a vanishing diffusive 
flux, d x s(x, t)\ x= i = 0, (Neumann boundary condition), 
i.e., a flat concentration profile. 

We assume now that a (non-reactive) tracer, e.g., 
some chemical species, flows diluted in the fluid which 
is injected into the porous column. Provided that the 
medium is sufficiently homogeneous, and that physical- 
chemical interactions of the transported species with 
the porous matrix and preferential flows can be ex- 
cluded EH > the solutes dynamics obeys an ADE with 
saturation-dependent coefficients 

d t sc(x, t) = -d x [u(s) - sD(s)d x ]c(x, t), (5) 

where c(x, t) is the solutes concentration, and s = s(x, t). 

The advection term u(s) is determined by the fluid 
flow, namely u(s) = j s (x, £), whereas the effective disper- 
sion coefficient D(s) accounts for the effects of mechanical 
dispersion and molecular diffusion mechanisms Q . Note 
that Eq. is linear, although knowledge of the satura- 
tion s(x,t) is required in order to determine c(x 7 t), i.e., 
problems ([5]) and ([4]) are inherently coupled. Owing to 
this coupling and to the nonlinearities, the flow patterns 
are not homogeneous, so that the evolution of the solutes 
dynamics usually displays non-Fickian features, such as 
long tails and non-Gaussian shapes. We will discuss this 
point in detail in the next Section. When the satura- 
tion level is uniform within the column, i.e., s(x,t) = sq, 
Eq. (J51) reduces to a standard ADE with constant coeffi- 



FIG. 4: Saturation level spatial profiles s(x,t) at fixed time 
t — 1 . Numerical integration is displayed as solid lines, Monte 
Carlo simulation as symbols: Ko = 0.4 in the left layer and 
fto = 0.05 in the right layer (crosses); kq = 0.05 in the right 
layer and ko = 0.4 in the left layer (circles). The other coef- 
ficients are: v = 0.2, vo = 0.3, and \ = 0.1. 

cients and contaminant transport becomes Fickian. 

Concerning initial and boundary conditions for the so- 
lute species, in the following we will assume that contam- 
inant release occurs within a given time interval to < t < 
t c , and that during this time span the pollutants concen- 
tration at the column inlet has a constant value Co, i.e., 
c(0,io < t < t c ) — co. Such finite-extension contaminant 
spills arc commonly encountered in environmental reme- 
diation problems [3J]. Before injection, there is no con- 
taminant within the column, i.e., c(x,t < to) = 0. The 
boundary condition for the concentration at the outlet is 
dictated by the outlet boundary for the flow, i.e., it must 
be a Neumann boundary condition, d x c(x, t)\ x= i = 0. 

III. A RANDOM WALK APPROACH 

Flow-transport equations ^ and © share a similar 
structure, and can be both written in conservative form 
as 

dt0p(x,t) = -d x [0q - 6dd x ]p{x,t), (6) 

for the evolution of the field p(x,t). The coefficients 
9 = 6(x,t), q — q(x,t), and d = d(x,t) are generally 
space and time dependent. This conceptual picture al- 
lows Eqs. (QJ and ([5]) to be expediently solved by re- 
sorting to a random walk formulation. The key idea is to 
think of the evolving (fluid or contaminant) plume, whose 
dynamics is described by Eq. ©, as being composed of a 
large number of particles performing stochastic trajecto- 
ries in the traversed porous medium. Then, the quantity 
p(x, i) > can be given a probabilistic interpretation 
(up to a normalization factor), i.e., p(x,t) represents the 
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FIG. 5: Breaktrough curve (circles) and saturation level 
(crosses) at the column outlet, normalized so to have unit 
maximum. Solid lines correspond to numerical integration, 
symbols to Monte Carlo simulation. A discrepancy between 
the profiles is apparent for non-constant velocities v. The 
coefficients are: ko = 0.05, v — 0.4, vq = 0.3, and \ — 0.1. 



probability density of finding a walker at a given position 
x, at time t. 

Consider an ensemble of N particles at positions Xj (t) , 
j = 1,...,N. Assume that the stochastic dynamics of 
each walker is governed by a Langevin-type equation 



Xj (t + r) = xj (t) + A(x, t)r + yj2B(x,t)r£j , (7) 

where A(x,t) is the drift coefficient, representing the av- 
erage velocity, B(x,t) is the diffusion coefficient, and r a 
(small) time step. The quantity £j is a white noise with 
zero mean and unit variance. The coefficients A(x,t) 
and B(x,t) completely define the properties of the mi- 
croscopic particles dynamics. Remark that Eq. |(7J) de- 
scribes a Markovian (memoryless) process: although the 
evolution of a single trajectory may depend on the others 
(i.e., the process can in general be nonlinear), knowledge 
of particles positions at time t is sufficient to determine 
the displacements at the following step t + r. It can be 
shown that the ensemble-averaged density P(x, t) of a 
particles plume obeying Eq. 0, i.e., 



P(x,t) = (£S[x-Xj(t)]), 



(8) 



satisfies the (in general nonlinear) Fokker-Planck equa- 
tion 

d t P(x, t) = -d x [A(x, t) - d x B(x, t)} P(x, t). (9) 

Nonlinearities arise when A = A{P) and/or B = B(P), 
i.e., when the coefficients depend on particles concentra- 
tion. 

Then, if we want to identify the walkers in Eq. (J7J) with 
the microscopic dynamics underlying Eq. ©, we must 



properly assign the drift A(x, t) and diffusion B(x,t) of 
the stochastic process. In other words, we must impose 
A(x, t) and B(x, t) so that knowledge of P(x, t) provides 
information on the quantity p(x,t): this in turn estab- 
lishes a link between the physical variables 9(x, t), q(x, t), 
d(x,t) and the parameters A(x,t) and B(x,t). 
Letting [U 



A(x,t) = q(x,t) + 



1 



6(x,t) 



d x 9{x,t)d(x,t) (10) 



and 



B(x,t) = d{x,t), 



it is easy to prove that we can identify 
p(x,t)=P(x,t)/6(x,t), 



(11) 



(12) 



up to a normalization factor, which provides the desired 
link. 

Let us address Eq. (j4} first. In this case, the particles 
that stochastically travel in the porous medium represent 
fluid parcels that progressively change the saturation dis- 
tribution in the traversed region [l2| . The nonlinearity of 
the governing equation arises from the fact that the ad- 
vection and diffusion coefficients depend both on s(x,t). 
Then, determining the evolution of the saturation profile 
at time t + r requires preliminarily knowing the satura- 
tion profile itself at time t. In terms of random walks, 
this implies that particles positions at the following time 
step can be updated once the positions of all the parti- 
cles at the current time have been determined. In other 
words, particles trajectories are correlated via the satu- 
ration, as the advection and diffusion coefficients depend 
on the fluid saturation, s. At each time step r, s(x,t) 
is first computed on the basis of particles positions at 
time t, then time is updated t = t + r and particles are 
displaced. From the previous considerations (the param- 
eter 9 is assumed to be constant and can be simplified), 
it follows that 



where 



and 



x °(t + T)=x s j (t)+ f x s + a s Zj, 
[v(s) +8 x k(s)]t 



(13) 
(14) 
(15) 



In the hydrodynamic limit r — > 0, the ensemble-averaged 
fluid parcels profiles P s (x,t) converge to the solution 
s(x,t) of Eq. 0}. 

Once s(x, t) has been obtained at each time step, the 
evolution of the concentration profile in Eq. §E§ can be 
determined from a second ensemble of particles repre- 
senting the pollutant parcels. This (linear) random walk 
must obey 



x c Jt + r)=x c Jt)+^ c + a^j, 



(16) 
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FIG. 6: Contaminant concentration profiles c(x, t)s(x, i) at a 
fixed time t = 1, for step injection from to — to t c — 0.1. 
The simulation parameters are ko = 0.05, \ — 0.4, vo = 0.3, 
v = 0.2, Do = 0.1, and <f> — 2. Monte Carlo simulations 
are shown as symbols (circles at t = 0.5, crosses at t = 1), 
numerical integration as solid lines. 



FIG. 7: Contaminant profile c(£,t)s(£,t) measured at the col- 
umn outlet, as a function of time, for step injection from 
to = to t c = 0.1. The simulation parameters are ko — 0.05, 
X = 0.4, v = 0.3, v = 0.2, Do = 0.1, and cf> = 2. Monte Carlo 
simulations are shown as symbols, numerical integration as 
solid line. 



where 



and 



u(s) 



-d x sD(s) 



y/2D( 8 y> 



(17) 



(18) 



Finally, we can identify c(x,t) = P c (x,t)/s(x,t), where 
P c (x,t) is the ensemble-averaged concentration of the 
contaminant walkers. 

In principle, the random walk schemes defined above 
can be used to determine fluid saturation and contami- 
nant concentration profiles for an arbitrary choice of the 
time and space dependent coefficients. In practice, how- 
ever, because of the sharp gradients and steep profiles 
resulting from the nonlinearities in Eq. special care 
is needed in the choice of the numerical values for the 
time step, r. Also, the evaluation of the space deriva- 
tives in the drift terms of the random walk is ill-defined 
for abrupt jumps (discontinuities) in the equations pa- 
rameters [2l|, [22|, [23| • In all such cases, it is convenient 
to resort to the ad hoc scheme originally proposed in [2l[ 
for particle transport in composite porous media. For 
instance, the random walk for the case of fluid parcels 
would read 



x ,■ 



„j(t + t) = x°{t) + v(s)t + a^x^t) + Ax%-, (19) 
where 



The expression for the case of contaminant particles 
(with non-constant 9) is more involved and reads 

x°(t + r)= x^t) + + -^U c [^(i) + A.t%, (21) 

s ys 



where 



and 



<7 C = y/2D(s)t 



A.T C = —<T C {x°(tM r 



(22) 



(23) 



Ax s = a s [x s JtMr 



(20) 



In [2l[ it has been shown that schemes (|19p and (p?Tj) 
are equivalent to (fT3|) and (fl"6|) . respectively, when coef- 
ficients are sufficiently smooth. 

Finally, just as for the governing equations above, 
the random walk schemes must be supplemented by the 
appropriate boundary and initial conditions. For each 
scheme separately, at time t = a given (large) number 
of particles is attributed to each dx along the discretized 
domain representing the column, so to reproduce the ini- 
tial saturation distribution and the contaminant concen- 
tration profile. A constant saturation level (or concentra- 
tion) at the inlet is imposed at each time step by keep- 
ing the number of particles located at x = equal to 
some reference value, i.e., by replacing the walkers that 
have either come back to the reservoir (x < 0) , or moved 
towards the interior of the column (x > 0). The Neu- 
mann boundary condition at the outlet is imposed by 
applying a reflection rule to the diffusive component of 
the displacement, while particles advected past the outlet 
during the same time step are removed from simulation. 
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Note that a Dirichlet (absorbing) boundary condition at 
the outlet would correspond to removing each particle 
from the ensemble upon touching the column end. 

Given an initial particles configuration, these are dis- 
placed at each time step according to the rules prescribed 
above. First, the necessary coefficients are computed at 
assigned s(x,t) profile (which is known on the basis of 
flow parcels positions). This allows displacing the walk- 
ers at position Xj during the time step r. Then, the 
s(x, t) profile is updated. Finally, the walkers at position 
Xj are displaced and their profile updated. Up to a nor- 
malization factor, the spatial profiles are determined by 
ensemble-averaging the walkers locations at a fixed time; 
the breakthrough curves at a given position are deter- 
mined by counting the net number of walkers crossing 
that location at each time step. Given the nonlineari- 
ties and the coupling between the two schemes, the time 
step t must be chosen sufficiently small to achieve con- 
vergence. Moreover, the number of simulated particles 
must be sufficiently large to attain a good accuracy in 
the estimated profiles. 



IV. NUMERICAL SIMULATIONS AND 
COMPARISONS 

We compare now the Monte Carlo simulations of the 
random walk schemes proposed in Sec. Mil to the nu- 
merical solution of the governing equations for s(x, t) 
and c(x,i). In principle, the random walk schemes are 
very general and can account for arbitrary functional 
forms (even discontinuous) of the various coefficients. 
In the following examples, we will focus on a power- 
law scaling, which is commonly encountered in the em- 
pirical constitutive laws, such as the Van Genuchten or 
Brooks and Corey laws, and can usually fit experimen- 
tal data 0, Q3, El El 111- In particular, for the vol- 
umetric fluid content we will assume re(s) = kqs x and 
v(s) = vos u , kq and vq being some constant reference val- 
ues. The exponents \ an d v are material parameters and 
depend on the details of the microgeometry. Moreover, 
for the solutes concentration evolution we also assume 
power-law scaling, D(s) = Dqs^ , where Dq is a constant 
diffusion coefficient and <f> is the scaling exponent. 

The physical quantities that we examine in the fol- 
lowing are the spatial profiles at a fixed time, which 
are helpful in estimating the average displacement and 
spread of the fluid flow and of the contaminant species, 
and the breakthrough curves at the outlet of the column, 
which allow assessing the distribution of the times needed 
to travel from the source to the measure point. In the 
following, we assume that the column has unit length, 
6=1. Furthermore, we assume a constant saturation 
level s(0,t) = 1 at the inlet, and let the fluid flow infil- 
trate the column under the combined action of capillarity 
and gravity. 

Fig. [1] shows the influence of the power-law scaling 
exponent for the capillary diffusivity, \. on the spatial 



profiles of the fluid saturation, s(x, t), at a fixed time 
t = 1 (x = 0.1, 0.5, 1 and 2). Saturation profiles become 
steeper as x increases. Fig. [2] shows the effects of the 
exponent of the velocity scaling law, v, on the spatial 
profiles s(x, t): again, the profiles become steeper as v 
increases, although the variation is milder than for the \ 
variation. In both cases, the agreement between Monte 
Carlo simulation and numerical integration of the corre- 
sponding equation is excellent. 

In the context of underground contaminant transport, 
abrupt spatial variations in the physical properties of the 
traversed media may commonly arise [2ll. l23j. These, in 
turn, give rise to sharply varying (i.e., possibly discontin- 
uous) transport coefficients and strongly affect particles 
trajectories. We address one such cases by considering a 
spatially discontinuous kq'. in particular, we assume that 
at a given interface between two layers kq has a sudden 
step variation, while being constant in each layer sepa- 
rately. This situation is usually referred to as a macro- 
scopic heterogeneity (the two layers are thought to be ho- 
mogeneous at the local scale) [2g|. All the other param- 
eters are constant across the interface. The schemes p^|) 
and (f2"Tj) would be suitable to deal also with more in- 
volved cases, e.g., multiple heterogeneities. Recent ex- 
perimental results [11 suggest that modeling transport 
through a sharp interface may require skewed flux cor- 
rections [HI . In this work, we do not address this issue 
and assume that the fluid flux is adequately described by 
Eq. ([2]), so that the random walk schemes above hold. 
In Figs. [3] and 2] we display the breakthrough curves 
j s (£,t) as a function of time and the spatial saturation 
profiles s(x, t) at a fixed time, respectively. We compare 
the curves for the case of fluid flow passing first through 
the layer at high kq = 0.4 and then through the layer 
at low kq = 0.05 with those where flow is in the oppo- 
site direction. For the boundary conditions considered 
here, fluid flow reaches the outlet earlier in the former 
case (Fig. [3]). Remark that for the chosen boundary and 
initial conditions the breakthrough curve reaches the sat- 
uration value j s (i,t — > oo) = vq. The saturation profiles 
along the column show sharp gradients at the interface, 
while preserving continuity (Fig. 3]). In both cases, the 
agreement between Monte Carlo simulation and numeri- 
cal integration of the corresponding equation is excellent. 

Note that for the case of non-constant velocity v, the 
breakthrough curve j s {£,t) — v{s)s\ x= i docs not coin- 
cide (up to a normalization constant) with the satura- 
tion s(i,t) measured at the outlet. This is immediately 
apparent from Fig. where we show Monte Carlo simu- 
lations and numerical integrations of the two curves (for 
the same parameters). This point might be relevant while 
applying inverse problem techniques to the estimate of 
model parameters on the basis of experimental data. 

Finally, in Figs. [5] and [7] we display the spatial pro- 
files and the outlet values of the contaminant concen- 
tration, respectively, corresponding to a finite-duration 
step injection. The computed quantity is c(x, t)s(x, t), 
i.e., the product of concentration and saturation. The 



FIG. 8: Unsaturated media with trapping processes: fluid 
flow profiles P at time t = 2.25 10~ 2 . For all curves, a = 0.6 
and A = 1. Crosses: ao — 4, a = 0.25, bo = 0.2, b — 0.5. 
Dots: ao = 1, a = 0.25, bo — 0.2, b = 0.5. Squares: ao = 0.1, 
a = 0.25, b = 0.5, b = 0.5. Circles: a = 5, a = 0.25, 
bo = 0.5, b = 0.5. The corresponding numerical integration 
curves are plotted as solid lines. 



spatial profiles are visibly skewed (Fig. [5]) , and this be- 
havior is reflected in the long tail of the concentration 
measured at the outlet of the column as a function of 
time (Fig. [7]). These features result from the coupling 
with the saturation and from the nonlincaritics involved 
in the flow-transport processes, and could possibly ex- 
plain the heavy-tailed breakthrough curves reported in 
the literature (see, e.g., 0) for nonsaturated homoge- 
neous porous media. 



V. FLUID FLOW THROUGH UNSATURATED 
HETEROGENEOUS MEDIA 

So far, we have focused our attention on the case of 
unsaturated homogeneous porous materials, i.e., materi- 
als that do not display any significant degree of disorder 
in the pore-space geometry. The homogeneity hypothe- 
sis is mirrored in the Markovian nature of the associated 
random walks: all spatial sites are statistically equiva- 
lent, and particles sojourns have the same duration r 
at each of them, so that trajectories have no memory 
of past positions. On the other hand, it is well-known 
that porous media are actually characterized by hetero- 
geneities at multiple scales, which ultimately affect fluid 
and contaminant particles displacements [3, [2j| . 

Consider a fluid flow in a complex (possibly fractal) 
porous microgeometry. In such a situation, the fluid 
parcels tend to flow in preferential channels [19| , so that 
the distribution of the sojourn times at each site is nec- 
essaril y n onuniform, as suggested by experimental evi- 
dence [30]. A detailed account of Richards' equation in- 
adequacy to explain a number of fluid flow experiments 
can be found in [l7], Hi], HH, HH and references therein. In 




FIG. 9: Unsaturated media with trapping processes: fluid 
flow profiles at time t = 2.25 10 -2 . For all curves, a — 0.4 
and A = 1. Triangles (P) and dots (P m ): A = a (P m ) a 
and B = b (P m ) b , with a = 2, a = 0.4, b = 0.4, b = 0.6. 
Squares (P) and crosses (P m ): A = a {P) a and B = b (P) b , 
with ao = 6, a = 0, bo — 0.2, b = 0.5. The corresponding 
numerical integration curves are plotted as solid lines. 



recent years, some extensions of the Richards equation 
have been proposed to take into account these phenom- 
ena [HI, [13, Hj|- In particular, it has been conjectured 
that the wetting front in infiltration processes through 
nonhomogeneous aggregated media remains immobile for 
long time periods, and the structural hierarchy of the 
soil structure is a primary cause of non-Fickian dynam- 
ics fPh , [35j . An expedient means of incorporating such ef- 
fects into the random walk described by Eq. (J7J is to allow 
for the possibility of trapping events of random duration 
after each displacement (of duration r). While several 
hypotheses can be made on the interplay between dis- 
placements and retention times, each corresponding to a 
distinct conceptual picture of the underlying physical sys- 
tem, here we follow the lines of a continuous time random 
walk a ppr oach called fractal Mobile-Immobile Model (f- 
MIM) 36J, which suitably generalizes the discrete-time 
process defined by Eq. ([7]). 

Within this framework, it is commonly assumed that 
trapping times between displacements obey 'fat-tailed' 
(power-law) distributions ip s (t) ~ < _1_Q , a > 0. If the 
decay is sufficiently slow, it is not possible to single out 
a dominant time scale (i.e., the mean of the pdf is not 
defined) , and particles can thus experience a large varia- 
tion of sojourn times, hence a broad spectrum of effective 
velocities at each spatial site [14 1. 

More precisely, assume 



(24) 



< a < 1, -0 being a pdf concentrated on with 



survival probability = ft ip(t')dt' = Xt- a /T(l - 
a) + JC(t), JC being integrable. Here A > is a scaling 
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factor that defines the strength of the trapping events. 
While in principle A could depend on x,t [37|, and even 
on the particles concentration, here for sake of simplicity 
we assume that it is constant. 

Consider the random walk defined by Eq. ([7]). At the 
end of each displacement, particles wait at the visited 
spatial site for a random time obeying the pdf ip s (t): 
these sojourn times can be very long as compared to the 
time scale r. Because of trapping events, the number 
of jumps performed by each walker in a given time span 
may greatly vary: this is an effective means of describ- 
ing heterogeneous media. Wc denote by P rn (x,t) and 
P l {x, t) the density of mobile and immobile fluid parcels, 
respectively, and by P = P m + P l , the total density. 

Let Xj j7l be the position of walker j just after the n- 
th step, which begins at time ij.„. This 'mobile' step 
is followed by the n-th 'immobile' period, with duration 
T^^ a W n: where the W n are independent random variables 
drawn from ijj (so that the pdf of T X l a W n is ip s )- Then 
we have 

x^ n+ i = x hn + At + V2Bt£„, (25) 

and 

tj,n+l = tj, n + T + T^Wn. (26) 

As a particular case, when A = the second equation re- 
duces to t_7,n+i = tj^+T and we recover the homogeneous 
random walk defined by Eq. {?])■ When A and B are uni- 
form and constant, in the hydrodynamic limit (r — > 0) 
the particles dynamics above defines a linear fractal-MIM 
model [36l . |38| . a generalization of the standard (linear) 
MIM model [39|]. In this case, it can be shown via sub- 
ordination that P satisfies an equation akin to Eq. ([5]), 
except that the left-hand side is replaced by \dt + Xd"]P, 
d" being a Caputo derivative of order a [36| (details are 
provided in Appendix . 

Unsaturated flow in porous media can be character- 
ized by droplets, slug-flows, Darcy's flow, or all of the 
above. Correspondingly, A and B may depend on P or 
P m . In highly unsaturated materials, we may conjecture 
that smaller pores, previously wetted and full of trapped 
fluid, modify the surface properties of larger pores where 
mobile fluid flows, so that A and B may depend also on 
P % . For sake of simplicity, we have neglected here other 
possibilities: for instance, additional nonlinearities would 
be introduced at strong solutes concentrations [40] . In all 
such cases, it is more convenient to derive the governing 
equation for the fluid parcels density by resorting to the 
relation between P m (x,t) and P l (x,t), as in [38| . 

We make use of the ancillary pdf f(x,t) for a walker 
of just being released from a trap at x, at time t, and 
denote by r(x,t) a possible source term. Particles that 
are mobile at time t were either released from a trap, or 
came from the source at time t — t', with < t' < r. 
We assume that the diffusive step of the displacement 
occurs at the end of the mobile period (see the discussion 
in Moreover, during the time interval [t — t',t], the 



motion of a walker is determined by the velocity field A: 
we denote the travelled distance by u x ,t,t'- 
Hence, we have 

P m (x,t)= f [f + r}(x-u x , t , t ,,t-t')dt', (27) 
Jo 

which implies 

[f + r]{x,u Xi t, T ,t-T) =r- 1 P m (x,t) + 0(r) (28) 

provided that A (hence u) and f + r are smooth [fnj . 

The density P l also depends on / + r: any immobile 
particle at time t was trapped at a time t' < t, at the 
end of a mobile step that began at time t — t' — r and 
involved a single diffusive step V2£>t£. Denoting by y 
the amplitude of this latter step, we have 

P*(M) - f yeR dyfo dt'^t'/r^^iy) 

[f + r](x-y- u x -y,t-t>,<r, t-t'- t). (29) 

Here ^(i'/r 1 /") represents the probability for a trap- 
ping duration to be larger than t', and f^rg^iy) = 

l/V2Brip(y/V2Br) denotes the pdf of £, ip being the 
normal distribution. Moreover, B may be nonuniform, 
and depend on the starting point of each jump (x — y in 
Eq. (f2"9"| ). In the following, we will assume that A and 
B depend on (x,t) either directly, or because they are 
functions of densities such as P m , P l or P. 

Denoting time convolutions (in R + ) by *, i.e., F * 
G(t) = f*F(t - t')G(t')dt', and recalling that * is 
bounded by 1, wc make use of approximation (|28p in 
Eq. ([29]) and obtain 

lyeR Pm ( x ~ V> ^^m^vV^dy + 0(r). (30) 
In the hydrodynamic limit, 

(due to Lemma 4 of Appendix |C|) . and the (time) convo- 
lution of kernel T~ 1, I'(t/T 1 / Q ) converges to the fractional 
integral A/g^_ Q (see Appendix IX)) . hence 

p\x 1 t) = \i 1 ; + a p m {x 1 t), (31) 

which provides a relation between P m and P l . We can 
identically rewrite as 

P(x, t) = [Id + \$- a ]P m {x, t). (32) 

Then, the inversion of the term Id + A/g yields [38| 

P m (x,t) = [Id + AJ 1 i + a ]- 1 P(a:,i). (33) 

The hydrodynamic limit of the walkers flux T(x, t) (cor- 
responding to the dynamics in Eqs. (f2"5)) and (|2"6"]) ) is de- 
rived in Appendix [B] and reads 

T=\A- d x B]P m , (34) 
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This follows from the fact that particles contribute to the 
flow only when being in the mobile phase P' m . Finally, 
mass conservation dtP = —d x T + r implies 

d t P = -d x [A-d x B]P m +r, (35) 

which is the desired governing equation for transport 
processes with trapping events [52||. Equation (f3"5)) . 
rather than being simply postulated as a phenomeno- 
logical 'fractional derivatives generalization' of the stan- 
dard Richards equation, has been here derived as the 
hydrodynamic limit of an underlying nonlinear and non- 
Markovian stochastic process with a definite physical 
meaning. Remark that Eq. (|35[) is in principle nonlin- 
ear, as the coefficients A and B may depend on P, P m 
or P 1 , as is the case for unsaturated flows. Moreover, 
Eq. (f35|) contains a memory kernel (via the relation (|32|) 
between P m and P): in other words, because of the 
slowly-decaying time kernel, knowledge of the past his- 
tory of the particle is required in order to determine its 
future displacement. Therefore, the fluid parcels density 
is affected at the same time by nonlinearities and mem- 
ory effects: the relative strength of these components ul- 
timately determines the fate of the flow in the traversed 
media. 

Observe moreover that Eq. (|35|) is written in 'Fokker- 
Planck form', i.e., with a term of the kind d x d x BP m . If 
we were to postulate a conservative (Fickian) formulation 
for the flux, corrections to the drift coefficient A should 
be introduced, in analogy with the case of the standard 
Richards equation for homogeneous media. 

The behavior of fluid flow with nonlinear coefficients 
and retention times is illustrated in Figs. [5] and [5J In par- 
ticular, we proceed to compare the Monte Carlo simula- 
tions of the random walks described by Eqs. ([26)1 and ([25]) 
(in the hydrodynamic limit) with the numerical integra- 
tion of the governing Eqs. (|3ip and (|35|) . Monte Carlo 
simulations proceed along the same lines as for homoge- 
neous media. Concerning the numerical integration, wc 
found more expedient to recast Eq. (|35[) in the equivalent 
formulation 

[d t + \D?]P m = —d x [A - d x B]P m + r, (36) 

-D" being a Riemann-Liouville fractional derivative, de- 
fined in Appendix |A"1 Once Eq. (|36p has been solved for 
P m , P is easily computed from Eq. (f32)) . We discretized 
Eq. (|3"6"]) according to a semi- implicit scheme as in [38j |. 
so to avoid possible instabilities connected with nonlin- 
earities. 

In Fig. [5] we display the total fluid flow density P at a 
given time, for a different values of the coefficients. We 
make the hypothesis that A and B depend on the mo- 
bile phase P m , with a power-law scaling A = ao(P m ) a 
and B = bo(P m ) b (similarly as done for the homogeneous 
media). The initial condition is a fluid pulse located at 
xq = 1/1. Absorbing boundary conditions are set at ei- 
ther end of the medium. In Fig. [5] wc display the total 
fluid flow density P as compared to the mobile density 



P m , when the parameters A and B separately depend 
(with a power-law scaling) on the mobile or total fluid 
density. Boundary and initial conditions are the same as 
in the previous example. Both figures show a very good 
agreement between Monte Carlo simulation and numeri- 
cal integration. 

Finally, replacing Eq. (|2"4")l with ip s (t) = Tipit/r) yields 
P % = XP m when ip has a finite average A, i.e., when the 
pdf decays sufficiently fast at long times. In this case, 
Eq. (35J) holds with P m = (1 + A) _1 P, and we have 

(1 + \)d t P = -d x [A - d x B]P + (1 + A)r, (37) 

which is a nonlinear Fokker-Planck equation with a re- 
tardation factor A [37l ] . 

VI. DISCUSSION AND CONCLUSIONS 

In this work we have addressed nonlinear coupled flow 
and transport processes in variably saturated porous me- 
dia. After briefly recalling the governing equations for 
homogeneous materials, we have shown how to describe 
fluid and solutes parcels trajectories by resorting to a uni- 
fied random walk framework. Monte Carlo simulations 
of the underlying microscopic particles dynamics have 
been successfully compared to the numerical solutions of 
the governing equations. The random walk approach has 
turned out to be a flexible tool, which allows dealing with 
nonlinear and/or discontinuous transport coefficients. 

Then, we have introduced a nonlinear f-MIM continu- 
ous time random walk scheme aimed at describing fluid 
flow through variably saturated heterogeneous media. 
The effects of spatial heterogeneities are mirrored in the 
possibility of long (power-law distributed) sojourn times 
of the flowing particles at each visited site. The cor- 
responding governing equations have been derived and 
their numerical integration has been then compared with 
the Monte Carlo simulations of the walkers dynamics. 
We remark that the proposed f-MIM generalization of 
nonlinear transport equations is not unique in any re- 
spect. Indeed, other possible approaches have been illus- 
trated, e.g., in fill. |42| and in [43| by means of a gener- 
alized Montroll- Weiss master equation with a jump pdf 
depending on the walkers density. 

Although focus has been given to fluid flow through 
heterogeneous unsaturated materials, the proposed non- 
linear f-MIM scheme is fairly general and can be straight- 
forwardly applied to the coupled solutes transport prob- 
lem, as well. Similarly as fluid parcels can be affected 
by the irregular geometry of the traversed material, the 
solutes concentration can also experience the effects of 
nonhomogeneities, due, e.g., to chemical-physical ex- 
changes of the solute species with the surrounding en- 
vironment [iijl . In analogy with the case of fluid flow 
dynamics, we may then take into account these contri- 
butions by introducing a power-law distribution ip c (t) ~ 
i -1- ' 3 , with (3 > 0, for the waiting times of contaminant 
particles between displacements. This pdf characterizes 
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the sorption times of the transported species within the 
medium. In general, there is no reason to suppose that 
the exponent f3 coincides with a. A similar distinction 
between flow- and solutes-induced retention times has 
been previously introduced in [3, 0, EH| within a (lin- 
ear) CTRW framework. 

In summary, the solutes concentration in heteroge- 
neous variably saturated materials may display devia- 
tions from standard Fickian behavior due to three con- 
current processes: i) the space- and time-varying satu- 
ration profile within the medium (nonlinear effects), ii) 
the spatial heterogeneities experienced by the fluid flow 
(memory effects), and Hi) the spatial heterogeneities ex- 
perienced by the solutes (memory effects) . As experimen- 
tal data such as contaminant profiles or breakthrough 
curves are usually limited and/or affected by measure- 
ment noise, distinguishing the effects of spatial hetero- 
geneities on flow and transport processes separately is an 
highly demanding task, and many research efforts have 
been recently devoted to this aim (see, e.g., Q and ref- 
erences therein). At present, it is therefore an open ques- 
tion whether these distinct contributions could be sepa- 
rately analyzed on the basis of measured data, or rather 
described by an effective CTRW model: we will discuss 
in detail this topic in a forthcoming paper. 

As a final remark, note that in this work we have 
confined our attention to the migration of nonreacting 
(passive) species. Nonetheless, all the random walk algo- 
rithms introduced here could be easily extended so as to 
describe the transport of radionuclides, by computin g th e 
decay time before simulating the particle trajectory [461 ] . 

APPENDIX A: FRACTIONAL INTEGRALS AND 
DERIVATIVES 

The fractional integral Iq + J of order a > is 

is, + m = Y^f\t-t) a - 1 m>&, 

which is a generalization of the usual multiple integrals 
to arbitrary (positive) order [13, EH . 

The Caputo fractional derivative d"f of order n < a < 
n + 1 is 

a?/(t) = i^-^fit), 

n being an integer [53, [5(| ■ The Riemann-Liouville 
fractional derivative Df is instead defined as 

Dff(t) = 9 t "+ 1 7 "+ 1 -«/(t) 

for n < a < n + 1. Note that the Riemann-Liouville 
derivative applies to slightly more general functions than 
Caputo's, and when both can be applied (i.e., for dif- 
ferentiable functions /) we have Dff{t) = dff{t) + 
t~ Q /(0+)/r(l - a) for < a < 1. 

In Section [V] we use the following Lemma. 



Lemma 1. (i) Suppose $ is a positive- valued, de- 
creasing function defined on R + , with ^(0) = 1 and 
= t- a /T(l - a) + JC{t), with K integrablc and 
< a < 1. Then, the convolution of kernel t" 1 \I'(</t 1 / q ) 
converges to Iq + in L P (R + ) 1 with 1 < p < +oo, when 
r -> 0. 

(ii) If \1/ is integrable, then the convolution of kernel 
t^^^t -1 ) converges to J R+ ^(t)dtld. 

Proof. Due to the above definitions, the convolu- 
tion of kernel iT(l - a)- 1 {t/T 1 / a )- a is exactly J*"". 

Moreover, the convolution of kernel T _1 / a /C(t/r 1 / a ) 
(as a mapping of L P (R + )) is an approximation to 
J R+ JC(t)dtld [43, |48|, so that the convolution of kernel 
t^/C^/t 1 /") tends to zero when r — > (due to a < 1), 
which proves point (ih Point (ii) is & direct consequence 
of the reference [13, HI] concerning approximations to Id. 

Remark. With <pe(y) = l/£ip(y/£), the space convolu- 
tion (denoted by*, with, f*g(x) = f R f(x — x')q(x')dx') 
of kernel ipi converges to Id when I — > pTl , [48| , in LP , 
and we have ipg * G(x) — > G{x) pointwise when G is a 
continuous function. 



APPENDIX B: PROBABILITY CURRENT 

We would like to prove Eq. (j34|) . Let us first introduce 
the probability current of the random walk defined by 
Eqs. and (fS5|) . The current is the probability for 

a walker to cross a given point x to the right during 
a time interval dt, minus the probability of crossing to 
the left, divided by dt. Upon multiplication by the total 
number of walkers involved in the random walk, it yields 
the average number of particles that cross x per unit 
time. 

Note that a walker that is not at the end of the current 
mobile period crosses x during [t — dt, t], provided that it 
was released from the trap at time t—t' (with < t' < t), 
between x — u x% t,t' ~ A(x — u x ^,t' 7 1 ~ t')dt and x — u x j,t' 1 
if A > 0. This occurs with probability 

la if + r]{x- u x , t< t>,t - t')A(x - u x , t , t >,t- t')dtdt' 
~ P m (x,t)A(x,t)dt 

when r — > 0. Moreover, a walker crosses 1 by a 
diffusive step with probability J y>a [f + r](x — y — 

u x ,t,T,t — r)$(y j \/2tB{x — y))dy per unit time, where 
$(z) = f+°° <p(y)dy. The quantity $(y/V2rI) repre- 
sents the probability for a given step to be larger than 
y. This latter probability approximates r _1 j y>0 P m (x — 

y, t)${y/ \J2tB{x - y))dy when r -> 0, due to Eq. (Un- 
collecting finally contributions of jumps to the left, the 
probability current of the random walk in Eqs. (|26|) 
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and (25]) is 

P m A(x, t) + r- 1 J y>Q P m (x - y, t)$(y/y/2TB(x-y) - 
P m {x + y, t)<5>{y/^2TB{x + y))dy. 

Appendix [C] shows that the above integral converges to 
—d x BP m J R y 2 ip(y)dy, due to the rapid decrease at in- 
finity of if and <J>. Hence, in the hydrodynamic limit the 
probability current is given by Eq. (|34"f . 

APPENDIX C: TECHNICAL LEMMA 

Lemma 2. Let $ be a diffcrcntiablc function, in- 
tegrable over R + , positive and decreasing. Suppose 
also that B, as a function of x, has a derivative 
B' satisfying \yB'{x + y)/B(x + y)\ < 2 every- 
where. Suppose then that G is an integrable func- 
tion whose derivative is uniformly bounded, and 
set J± = t- 1 $ y>Q G(x ± y)$(y/y/2TB(x±y))dy. 
Then, the quantity —1+ + I— converges to 
-4d x (G(x)B(x)) / + °° z<S>(z)dz when r -> 0. 

Remark. We have —4 J + °° z$(z)ciz = f R z 2 ip(z)dz, 
given that (p(z) = — $'(z). 



and 



Proof of Lemma 2. We first show that (i) — J+ + J_ is 
:z)<I>(z)(iz, i* 1 being a de 
0. Then, («i) we apply Lemma 



of the form 2s" 1 J + °° F(ez)$(z)dz, F being a derivable 



function satisfying F(0) 
3 further below. 

In view of (i), let us fix x. The hypotheses of the 
Lemma allow for a change of variables z = g±{y)/V2r 
in J-j-, with g±(y) = y / \J B(x ± y). The inverse of g± is 
/i±, with 



<4(y) = l/y/B{x ±y)[lT yB'(x ± y)/(2B(x ± y))], 

fe±(Z)B , (a:±fe ± (Z)) 1 -l 
2B(x±/i ± (Z)) J ' 



and ft.±(0) = 0. With these notations, set e = \/2t and 

F{Z) = G(x-h-(Z))h'_(Z)-G(x+h + (Z))h' + (Z) (CI) 

so that letting y = h±(ez) in L± yields —I + + /_ = 
2s' 1 f+°° F(ez)$(z)dz. 

Now, to apply Lemma 3, note that F'(Q) = 
-2G'(x)B(x)+G(x)(h'!_(0) - ^'j(O)), and that h'l(Z) = 
a± ± b±, with 



a± 



± h' ± (Z)B(x ± h ± {Z)) h ± (Z)B'(x ± h±(Z)) 1 _ 1 



y/B( X ±h ± (Z)) 



with 



m h ± (z)B'(x±h ± (z)) ]2 2B(x ± ft±(Z)) 2 ' 



c± = h' ± (Z)BB'{x ± /i±(Z)) ± h±h' ± (Z)BB"{x ± /i±(Z)) 
T fc±(Z)h / ± (Z)fi /a (x±h±(Z)), 

which implies /4(0) = ±S'(x), so that F'(0) = 
—2(GB)'(x). Hence, Lemma 2 is a consequence of the 
Lemma 3 below, which itself is included in Lemma 3 

of [33. 

Lemma 3. Let $ be a differentiablc function, integrable 
over R + , and bounded. Then, for any integrable function 
F satisfying -F(O) = and whose derivative is uniformly 
bounded, the expression 



F(ez)$(z)dz 



2y/B(x±h±(Z)) 



2B(x±h±(Z)) 



converges to F'(0) / + °° z${z)dz when e -> 0. 

In Section [V] we use the following statement. 

Lemma 4- Let G be a continuous bounded 
function, with (p a pdf. Suppose also that 
B is as in Lemma 2. Then, / = J R G(x — 
y)<p(y/y/2TB(x - y))dy/yj2 T B{x - y) when 
r 0. 

As in the proof of Lemma 2, the change of variables 
g(y) = yj B(x — y) has an inverse, which we denote by 
h. Thus, we have I = J B G(x — h{\f2~Tz))ip(z)dz. Hence, 
the remark of Appendix IA1 proves the statement. 
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